統計勉強メモ2

お勉強メモ2 集団の比較について

In [1]:
import sys

import matplotlib.pyplot as plt
from matplotlib.pyplot import *
import numpy as np
from numpy import *
import numpy.matlib
from scipy import stats as sps

1. 集団の母平均比較

2つの集団の標本平均を比較して，本当に差があるのかどうかを調べたい．
集団に対して分かっている情報によって方法が変わる．
このページが体系的にまとまっているので引用
http://www.geisya.or.jp/~mwm48961/statistics/mobile/bunsan1_m.htm

母分散が等しい 母分散が異なる
データ対応が有る A A
データ対応が無い B C

A,B,Cのそれぞれのパターンについて比較方法が異なる．

1.1 データに対応がある場合

データの対応

という言葉の定義ははっきりしたものは見つけられなかった．

なので自分なりに定義

ある母集団$X$の部分集団$\tilde{X}$に対して，何らかのテスト$f,g$を行って得られるテスト結果$f(\tilde{X}),g(\tilde{X})$があるとき，テスト結果同士の関係性について「部分集団$\tilde{X}$による対応がある」と定義する．

この定義を用いて，今から考えたい問題は以下のように定義できる．

母集団$X$にテストを行った場合に想定される結果$f(X),g(X)$の差の調査を，部分集団$\tilde{X}$による対応があるテスト結果$f(\tilde{X}),g(\tilde{X})$を用いて行う．

方針

知りたいのはテスト$f,g$の差なので以下の確率変数$z$の期待値$E[z]$とその推定精度 $$z=E\left[f(X)\right]-E\left[g(X)\right]$$

ここで，新たな確率変数$y$を以下のように定義する． $$y_i \equiv f(\tilde{x}_i)-g(\tilde{x}_i)$$

このとき，確率変数$y$の期待値は知りたい確率変数$z$の期待値に一致するはず(要検証) $$E[y]=E[z]=\bar{y}$$

ここで気になるのは，確率変数$z$の期待値$E[z]$の推定精度=標本平均$\bar{y}$がどれくらいの推定精度か？つまり，$\bar{y}$の確率分布はどのようになっているか？
その答えはt分布．
前回のStaticsノートで書いたように，母分散が未知の場合の標本平均$\bar{y}$の確率分布はt分布を使えば求めることができる．

その結果，「95%の確率でテスト結果$f(X),g(X)$の差は[標本平均]±[推定精度]である」と主張できる．

実例
この問題を考える． http://kogolab.chillout.jp/elearn/hamburger/chap5/sec3.html

この問題では対応のある2つのテスト結果が与えられる．
目的は2つのテストの点数に差があるかどうかを調査すること．
そのテストに関して，以下の情報が分かっている．

  • 標本数$n$=8
  • 点数差の標本平均$\bar{y}=-4.38$
  • 点数差の不偏分散$u^2=17.41$

テストの点数差の期待値は標本平均$\bar{y}=-4.38$に等しい．
あとは，この点数差がどれだけの精度か？を調査する．

今回のケースでは点数差の標本平均$\bar{y}$は以下の分布に従う．

  • 自由度7のt分布をscale倍にスケーリングし．locationだけオフセットした分布
  • scale : $\sqrt{\frac{u^2}{n}}=1.48$
  • location : $\mu$

自由度7のt分布に関して，全体の95%が含まれるようなtの値は，t分布表より2.365と与えられる．
なので，点数差の標本平均$\bar{y}$が95%含まれるような境界の値は$\mu\pm(2.365*1.48)=\mu\pm3.5$で与えられる．
よって点数差の母平均$\mu$は95%の確率で$-4.38\pm3.5\iff-7.88\sim-0.88$の範囲に含まれる．
つまり，「このテストの点数差について95%の確率で少なくとも-0.88点の差はある」ことがいえる．

1.2 データに対応がない場合

知りたいのはテスト$f,g$の差なので以下の確率変数$z$の期待値$E[z]$とその推定精度 $$z=E\left[f(X)\right]-E\left[g(X)\right]$$

データに対応がない場合，$y_i \equiv f(\tilde{x}_i)-g(\tilde{x}_i)$の傾向を見ても意味がない(というかそもそもそんなデータ作れない) なので，以下の統計量を考える．

$$\bar{y}=E\left[f(X)\right]-E\left[g(X)\right]\\ =\bar{x_a}-\bar{x_b}$$

ここで，$x_a$,$x_b$はそれぞれテスト$f,g$を実施した集団の要素である．
このとき，確率変数$\bar{y}$の期待値は知りたい確率変数$z$の期待値に一致する(要検証) $$E[\bar{y}]=E[z]$$

ここでも気になるのは，$\bar{y}$の確率分布である．
どうやら標本平均の差は以下の分布に従うらしい

  • 自由度dofのt分布をscale倍にスケーリングし．locationだけオフセットした分布
  • 2つの集団の母分散が等しいとみなせない場合(Welchの近似)$$ dof=\frac{\left( \frac{u_a^2}{n_a}+\frac{u_b^2}{n_b} \right)^2}{\frac{u_a^4}{n_a^2(n_a-1)} + \frac{u_b^2}{n_b^2(n_b-1)}}$$
  • 2つの集団の母分散が等しいとみなせる場合$$dof=n_a+n_b-2$$
  • $scale=\sqrt{\frac{u_a^2}{n_a}+\frac{u_b^2}{n_b}}$
  • $location=\mu$

ここで，標本平均の差の標準偏差については分散の加法性から求めている．

分散の加法性 $$V[x+y]=V[x]+V[y]+Cov[x,y]$$

ここでは，$x=\bar{x_a}$,$y=\bar{x_b}$である． また，$\bar{x_a}$と$\bar{x_b}$は独立と仮定すれば$Cov[\bar{x_a},\bar{x_b}]=0$である．

Welchのt検定を使えば，2つの集団の母分散について何も知らなくても母平均の比較ができる． なので，母分散が等しいとみなせる場合の自由度の式は不要に思える． しかし，「統計的方法演習(草場郁郎，日科技連，1974，p.24)」では以下の記述がある．

実用上の問題としては分散に違いがあるときに平均値の差を論じても意味のないことが多い．まず，分散にさのある原因を追求し，処置をとって，等分散となってから平均値の差について考えるほうが良いことが多い．
Welchの方法を用いれば，この解答のようにF検定をしてからt検定をするという二重の手間をかけず，ただちに平均値の検定を行えるわけである．しかし，実際問題においては，分散の違いがあるかどうかという情報をF検定で得ることのほうが検定の手間をはぶくより意味のあることが普通である．

とりあえずWelchのt検定使っておこうみたいな安易な考え方は危ないかもしれない．

2．集団の母分散比較

母平均の比較と同様に，不偏分散の差に関しても何らかの分布に従うのではないか
→不偏分散の比が従う分布がある．その分布をF分布という．

不偏分散の比$\frac{u^2_a}{u^2_b}$は以下の分布に従う．

  • 自由度[dof1,dof2]のF分布をScale倍にスケーリングした分布
  • dof1=$n_a-1 \equiv \phi_a$
  • dof2=$n_b-1 \equiv \phi_b$
  • scale=$\frac{\sigma^2_a}{\sigma^2_b}$

また，不偏分散の比について，以下が成り立つ

  • 最頻値は以下の式で与えられる $$\frac{\phi_b (\phi_a-2)}{\phi_a(\phi_b+2)} \cdot scale$$
  • 平均値は以下の式で与えられる $$\frac{\phi_b}{\phi_b-2} \cdot scale$$

https://staff.aist.go.jp/t.ihara/f.html

http://bio-info.biz/statistics/distribution_f_distribution.html

In [2]:
na,nb=8,15
fa,fb=na-1,nb-1
sig2a,sig2b=2,9
scale=sig2a/sig2b
mua,mub=-1,3
m=5000

xa=np.random.randn(na,m)*sqrt(sig2a)+mua
xb=np.random.randn(nb,m)*sqrt(sig2b)+mub
u2a=np.var(xa,axis=0,ddof=1)
u2b=np.var(xb,axis=0,ddof=1)

u2r=u2a/u2b

# fdistribution
t=np.linspace(0,3,500)
f=sps.f.pdf(t,dfn=fa,dfd=fb,scale=scale)

# center of f
c=(fb*(fa-2))/(fa*(fb+2)) * scale
cq=np.median(u2r)
# mean of f
ave=fb/(fb-2)* scale
aveq=np.mean(u2r)

figure()
g=plt.hist(u2r,bins=100)
f=f/max(f)*max(g[0])
plot(t,f,'r')
plot([c,c],[0,max(g[0])],'c',label='mode')
plot([ave,ave],[0,max(g[0])],'y',label='mean (theory)')
plot([aveq,aveq],[0,max(g[0])],'m--',label='mean (experim)')
legend()
plt.show()

応用

不偏分散の比$\frac{u^2_a}{u^2_b}$を，自由度[dof1,dof2]のF分布の上2.5%,下2.5%の値で不偏分散比を除したものが，母分散比の95%信頼区間になる．

つまり，この区間内に95%の確率で母分散比が含まれるということ．もしこの範囲に1が含まれていたら，母分散の比に差があるとは言えない．

In [3]:
Fu,Fl=sps.f.ppf(0.975,fa,fb),sps.f.ppf(0.025,fa,fb)
rangetable=(u2r/Fu,u2r/Fl)
tf=np.logical_and(rangetable[0]< scale , scale <rangetable[1])
ptable=[sum(tf),m-sum(tf)]
figure()
plt.pie(ptable,autopct="%1.2f%%",colors=["gold","lightpink"] \
        ,labels=["in","out"],startangle=90,counterclock=False)
axis('equal')
plt.show()

3. 3つ以上の集団の比較

t検定を全ての組み合わせについてやるのは良くないらしい．

大きく分類すると一元配置(一因子)，二元配置(2因子)，多元配置(多くの因子) になる．

比較の結果分かるのは比較してる集団のどれかに有意な差がある，ということだけで，実際にどの集団同士に差があるのか？については多重比較という方法を使わないといけない．
根が深そうなので必要になったらその時勉強する．